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Abstract 

We discuss a new fifth-order, semi-discrete, central- upwind scheme 
for solving one-dimensional systems of conservation laws. This scheme 
combines a fifth-order WENO reconstruction, a semi-discrete central- 
upwind numerical flux, and a strong stability preserving Runge- 
Kutta method. We test our method with various examples, and 
give particular attention to the evolution of the total variation of the 
approximations. 
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1 Introduction 

In this paper we present a fifth-order, essentially non-oscillatory central- 
upwind scheme that is designed to solve systems of conservation laws of 
the form 

* + /(«)* = 0 . ( 1 . 1 ) 

Here q E is a p-dimensional solution vector and / is a p-dimensional flux 
function. The solution of (1.1) may become singular in finite time, which in 
turn requires a careful study when dealing with numerical approximations. 
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One approach to approximating solutions of (1.1) is to use high-order, 
non-oscillatory central methods, which were introduced in [14]. Central 
methods avoid approximating the solution of (1.1) at singularities of the 
solution, and so do not require solving Riemann problems. The resulting 
simplicity makes central schemes well-suited for systems and multiple di- 
mensions. Central-upwind schemes, introduced in [7] and refined in [5], are 
semi-discrete variants of central methods which have improved efficiency 
and less dissipation than fully-discrete central methods. Our work uses the 
numerical flux of [5], which we refer to as the KNP flux. 

In this work we combine the KNP flux with the fifth-order weighted 
essentially non-oscillatory (WENO) reconstruction of [3], and the five-stage 
fourth-order strong stability preserving (SSP) Runge-Kutta method of [17], 
which is based on [1]. This is the first time these particular ingredients are 
combined into one scheme. 

Fourth-order fully-discrete central schemes based on WENO reconstruc- 
tions were presented in [9, 12]. The total variation behavior of these 
methods was examined in [11], where numerical experiments indicate that 
though the WENO-based methods axe not total variation diminishing (TVD), 
they are total variation bounded. Fifth- and ninth-order fully-discrete cen- 
tral schemes are discussed in [15]. Third-order extensions of the KNP 
scheme can be found in [4, 6]. 

In this paper we investigate the evolution in time of the total variation 
(TV) of our scheme. The total variation is defined for a discrete solution 
Uj as TV (u) := J2j l^j+i — u j\- In the case °f systems TV is defined as 
the sum of the TV over the components. A scheme is called TV bounded 
(TVB) in 0 < t < T if TV (u) < K for fixed K > 0 which depends only on 
initial conditions, and Vn and At such that nAt < T. In the scalar case, 
if a scheme is TVB then there exists a convergent subsequence in L\ oc to a 
weak solution of (1.1), which turns into strong convergence if an additional 
entropy condition is satisfied (see [8]). Our numerical experiments indicate 
that our method is TVB, providing evidence of the convergence of the 
method. 

The structure of this paper is as follows: In Section 2 we present our 
fifth-order central-upwind scheme, - summarizing the derivation of the KNP 
flux in Section 2.1. The WENO reconstruction is presented in Section 2.2 
and the fourth-order SSP-RK method in Section 2.3. Section 3 presents 
the results of a number of numerical tests of our method. We test both the 
accuracy and the evolution of the total variation of the resulting approxi- 
mations. 
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2 The Numerical Scheme 

We briefly summarize the components we use to construct our fifth-order 
central-upwind scheme; the numerical flux from [5], the reconstruction from 
[3] and the ODE solver from [17]. This is the first time these particular 
ingredients are combined into one scheme. 


2.1 The KNP Flux 


Throughout this section we assume a one-dimensional grid {xj} with con- 
stant spacing Ax. We define x J± i := Xj±| Ax and the cell Ij — jx^.i, x^+i 
For any function / (x) we use the notation fj := f ( Xj ). The cell average 

of q in the cell Ij is given by qj := f x q (x) dx. 

I 

We assume that the cell-averages q^ are known at time t n . The first 
step in the derivation of the approximate solution is to generate a piecewise- 
polynomial reconstruction from these cell- averages. Such a global recon- 
struction is defined as q (x) = qj (x) Xlj ( x ) > where Xi, (#) is the char- 
acteristic function of Ij , and qj (x) are polynomials of a suitable degree. 

In each cell Ij the reconstruction qj(x) should be conservative, i.e. 
~Kx Ix j+ ^ 0*0 dx ~ Qji sth-order accurate, so qj (x) — q(x) -f 0(Ax s ) 

for x € Ij, and non-oscillator y. Given such a reconstruction, we denote the 
point- values of q at the interfaces of the cell Ij by q+ + \ := qj + 1 ^x J+ i j 

and q~ +i -qj (*,+*). 

The local speeds of propagation of information from the discontinuities 
at the cell interfaces, , are given by 





min 




Here, A j denote the eigenvalues of the Jacobian of / evaluated at 

x j+ 1 . These local speeds of propagation axe then used to determine evolu- 
tion points that are away from the propagating discontinuities. An exact 
evolution of the reconstruction at these evolution points is followed by 
an intermediate piecewise polynomial reconstruction and finally projected 
back onto the original cells, providing the cell- aver ages at the next time- 
step <7™ +i . Further details can be found in [5] . _ A semi-discrete scheme is 
obtained in the limit as At — ► 0, yielding the KNP central- upwind scheme 


dqj _ jVl 

dt x j+\ ~ x j- i 


( 2 . 1 ) 
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The numerical flux in (2.1) is given by 



The accuracy of this scheme is determined by the accuracy of the recon- 
structions and the ODE solver. 

It is straightforward to generalize this scheme to higher dimensions, us- 
ing dimension-by-dimension reconstructions. Care must be taken, however, 
to use higher-order quadratures in the derivation of the KNP flux in higher 
dimensions to maintain accuracy. See [6] for a third-order example. 


2.2 The Fifth-Order WENO Reconstruction 

Weighted, essentially non-oscillatory (WENO) reconstructions [3, 13] are 
based on the essentially non-oscillatory (ENO) reconstructions of [2, 16]. 
ENO schemes choose the stencil that provide the least oscillatory recon- 
struction. WENO schemes weight all stencils so that accuracy is gained in 
smooth regions while avoiding crossing discontinuities. 

We use the fifth-order WENO reconstruction of the point-value 
given in [3]. We begin with the three quadratic reconstructions on three- 
point stencils 

2 

i ~ ^ ^ a rQj+k+r- 2? (2-2) 

r=0 

where k ranges from 0 to 2 and the coefficients a*, given in Table 2.1, 
are defined so that (2.2) approximates q 1 ^ with third-order accuracy. 
The WENO reconstruction is then defined as the convex combination 




(e+s>y 


The constants C k = {1/10,6/10,3/10} are defined so that Y^k=o C k Qj + i 
approximates q (xj + i ^ with fifth-order accuracy. S k is a smoothness mea- 
sure which is large when q k ' i has large variation. S k approximates the 
1 / 2 -norm of the first two derivatives of g, and is given by 

Sj ~ Y 2 ~ 2 “ 2gj-i + Qj) + — {qj-2 — 4 Qj—i T 3 qj) 

13 1 

s } = ^2 ^ j ~ 1 ” + 5 ?+ 1 ) 2 + 4 fe- 1 “ 4 ?+ i) 2 

Sj = ^ (5j “ 2g^ + i -f gj+ 2) 2 + - (3 qj - 4g J+ i + gj+ 2) 2 • 


(2.4) 
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Following [3] we take e = 10 6 . The reconstruction of q 3 _ i on the stencil 
centered at Xj is obtained by symmetry. For more details consult [3], 


k 

r = 0 

r = 1 

r = 2 

0 

1/3 

-7/6 

11/6 

1 

-1/6 

5/6 

1/3 

2 

1/3 

5/6 

-1/6 


Table 2.1: Coefficients a* for the quadratic reconstructions (2.2). 


2.3 The Fourth-Order SSP Runge-Kutta ODE Solver 

We use an optimal, strong stability preserving fourth-order accurate five- 
stage Runge-Kutta solver from [17]. We write this method as an explicit 
5-stage Runge-Kutta solver for the ODE dq/dt = F (q) as 


9 (0) = 


s — 1 


Q {S) = +At ^F(q (l) )) , 

(2.5) 

i = 0 V 


,"+1 = g (4) t 



where s = 0, . . . , 4 and 

a 0 ,o = 1-0, a li0 = 0.44437049406734, <*i,i = 0.55562950593266, 

c* 2 ,o = 0.62010185138540, a 2 ,i = 0, £* 2 , 2 = 0.37989814861460, 

a 3 ,o = 0.62010185138540, a 3 ,i = a 3 , 2 = 0, a 3 , 3 = 0.82192004589227, 

a 4 ,o = 0.00683325884039, a 4 ,i = 0, a 4 , 2 = 0.51723167208978, 

a 4;3 = 0.12759831133288, <* 4 , 4 = 0.34833675773694, 

0o,o = 0.39175222700392, 

/?i,o = 0, 0i,! = 0.36841059262959, 

02.0 = 02,i = 0, 02,2 = 0.25189177424738, 

03.0 = 03,1 = 03,2 = o, 03,3 = 0.54497475021237, 

04.0 = 04,1 = 04,2 = 0, 04,3 = 0.08460416338212, 

04,3 = 0.22600748319395. 

In our examples we use the CFL condition At = 0.9 min, where 
here A j denote the eigenvalues of the Jacobian of / evaluated at x 3 . 
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3 Numerical Examples 

We test our method on various examples, measuring both the accuracy 
and the evolution of the TV over time. When an exact solution is not 
available, we use a high-resolution reference solution computed using the 
second order KNP method as presented in [5], which is known to be TVD. 
Unless otherwise stated, this reference solution uses N = 5000 nodes. 

We first test our method with the scalar advection problem u t + u x = 0, 
u(x,t = 0) = sin 4 (7rx) on the periodic domain [—1,1] to T = 2. The 
relative L 1 - and L°°- norms of the errors are shown in Table 3.1. We also 
test accuracy with the Burgers equation = 0, u (x, t — 0) = 

3+sin (x) on the periodic domain [0, l7r] at T — 0.5, before shock formation, 
and at T = 2.5 after shock formation. The results are also shown in 
Table 3.1. Figure 3.1 shows the result at T = 2.5, as well as the change 
in TV over time for various resolutions. For this example the exact TV 
equals to 4 before the singularity formation at T = 1. 


N 

relative L 1 -error 

L 1 -order 

relative L°° -error 

L°°-order 

Linear Advection of sin 4 (7rx), T — 2 || 

100 

8.56xl0~ 4 

- 

2.73X10 -5 

- 

200 

2.58xl0 -i> 

5.05 

5.53 xKT v 

5.62 


6.68X10- 7 

5.27 

7.05X10- 9 



1.90xl0“ a 

5.14 

6.84X10- 11 



Burgers equation before the shock, T — 0.5 ] 

100 

1.09x10-° 

- 

1.08xl0- v 

- 

IKiCTI 

4.30x10-° 

4.66 

1.94X10" 9 

5.80 

400 

1.46X10 -9 ' 

4.88 

3.62X10- 11 

5.74 

800 

1.07xl0 _lu 

3.78 

8.38 xlO- 13 

5.43 


Burgers equation after the shock, T — 2.5 

■firm 


- 

1.74xl0 _s 

- 


HBsEESSU 

0.70 

1.40xl0- s 

0.31 


■ 

1.28 

4.98X10- 4 

1.49 


4.16X10- 4 

0.93 

2.61 xl0“ 4 

0.93 


Table 3.1: Relative L l and L°° errors for the advection equation and the 
Burgers equation 


Our next example is Burgers equation on the same domain with initial 
data u (x, t = 0) = 2 — sin (x) -f sin (2x). This example develops two shocks 
which eventually merge. Figure 3.2 shows the solution at T = 1.2 as well 
as the change in TV over time for various resolutions. 
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Burgers equation at T-2.5, N«1 00 Total variation for Burgers equation to T-2.5 



Figure 3.1: Results for the Burgers equation using the central- upwind 
scheme (2.1), (2.3) and (2.5). Left: the solution after shock formation 
at T = 2.5, exact solution, “o”: approximation. Right : the change in 
the TV of the approximation for N = 100, 200, 400, 800 nodes (from left to 
right) compared with the TV of a reference solution (the upper curve). 


2-shock Burgers problem at T-1 .2, N-1 00 Total variation for 2-shock Burgers problem to T-5 




Figure 3.2: Results for Burgers equation with initial data that develops a 
double shock using the central- upwind scheme (2.1), (2.3) and (2.5). Left : 
the solution after shock formation at T — 1.2, exact solution, “o”: 
approximation. Right: the change in the TV of the approximation for 
various resolutions compared with the TV of a reference solution. 
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Turning to systems, we consider the Euler equations 

( pu) +( fnf+p ) = 0, (3.1) 

\ E ) t \(E + p)uJ x 

with equation of state p = (7 — 1) (E — \ pu 2 ) and 7 = 1.4. We first apply 
our method to the Sod problem on the domain [0, 1] with initial data 


(p, u, E) 


(1.0, 0.0, 2.5), x < 0.5, 
(0.125,0.0,0.25), x > 0.5. 


(3.2) 


The results at T - 0.16 with N = 100 grid nodes is shown in Figure 3.3 
where we see that the shocks and the contact discontinuity are well captured 
by this scheme even at this low resolution. Figure 3.3 also shows the TV 
behavior of the approximation, compared with a reference solution. We 
see that the TV of the approximate solutions is greater than that of the 
reference solution, but damps with time suggesting that the approximations 
are TV bounded. 


Sod problem density at T-0.16, N-100 Total variation lor Sod equation to T*G.16 




Figure 3.3: Results for the Sod problem (3.2) using the central-upwind 
scheme (2.1), (2.3) and (2.5). Left: Density. reference solution, “o”: 

approximation with N = 100 nodes. Right : the change in the TV of the 
approximation for various resolutions compared with the TV of a reference 
solution. 

We now apply our method to the Lax problem on the domain [0, 1] with 
initial data 

(p, u, E) = 

The results at T = 0.16 with N = 100 and N — 400 grid nodes is shown in 
Figure 3.4. While the shock and the contact discontinuity are well-captured 


(0.445,0.311,8.928), x < 0.5, 
(0.5,0.0,1.4275), x > 0.5. 
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at low resolution, there is significant oscillation between the contact discon- 
tinuity and the shock for N = 100. Figure 3.4 also shows the TV behavior 
of the approximation, compared with a reference solution. Similarly to the 
Sod problem, we see that the TV of the approximate solutions is initially 
greater than that of the reference solution, but converges to the TV of the 
reference solution. 



Figure 3.4: Results for the Lax problem (3.3) using the central- upwind 
scheme (2.1), (2.3) and (2.5). Left: Density reference solution, “o”: 
approximation with N = 100 nodes, approximation with N = 400 

nodes. Right: the change in the TV of the approximation for various 
resolutions compared with the TV of a reference solution. 


For our final example we apply our method to the problem of a mach 3 
wave interacting with an acoustic shock on the domain [0, 1], see [16]. The 
initial conditions for this problem are 


f (3.857143,2.629369,10.3333), x < 0.1, 
\ (1 4- 0.2 sin (50x) , 0, 1) , x > 0.1. 


(3.4) 


The results are shown in Figure 3.5, compared with a reference solution 
using 20,000 nodes. We see the high resolution of our method and indica- 
tions that the TV of our approximations converges to that of the reference 
solution. 
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Acoustic problem density at T-0.18, N=400 Total variation tor acoustic problem to T=5 




Figure 3.5: Results for the acoustic-shock problem (3.4), showing the ap- 
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Approximate solution. reference solution, “o”: approximation. Right: 
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